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Abstract 



^t- : 
o : 
o ■ 

£N| ■ We introduce a numerical method for generating the approximating polynomials used in fermionic calcu- 

» I ' lations with smeared link actions. We investigate the stability of the algorithm and determine the optimal 

£1+ weight function and the optimal type of discretization. The achievable order of polynomial approximation 

' reaches several thousands allowing fermionic calculations using the Hypercubic Smeared Link action even 

with physical quark masses. 

0\ : 

^ ; 1 Introduction 

' The usage of smeared or fat links improves flavor symmetry for staggered fermions ^ [21 E] • Since the smearing 
f^i i contains a projection onto SU (3), the Hypercubic Smeared Link (HYP) action is not bilinear in the original thin 
link variables. Therefore, the explicit form of the fermion force is rather complicated making the Hybrid Monte- 
Carlo (HMC) and other molecular dynamics based algorithms with the HYP action practically unusable 1 . An 
update method based on a stochastic estimator [7J |S] can avoid this problem. The algorithm using improved 
stochastic estimators [S] requires polynomial approximation of functions of type x~~ a e p ( x >, where a > and 
p(x) is a low order polynomial. When the calculations are made at the small physical quark masses the order of 
these polynomials have to be in the range of the thousands. We introduce a numerical method to generate these 
high order polynomials and investigate the stability, optimal weight function and optimal type of discretization 
for the algorithm. In contrast to exact methods our procedure is very fast, stable up to thousands of orders and 
can be applied practically to functions of any type. 

The rest of the paper is organized as follows. In the next section we summerize the properties of the HYP 
action. In Sect.|3we describe the process of generating the approximating polynomials. We conclude in Sect.0] 



S ; 2 The HYP action 

The construction of the Hypercubic Smeared Link (HYP) action goes as follows First, the original thin 

links Ui^ are used to construct the set of decorated fat links, Vi >Kvp with a modified projected APE blocking 
step 



(1) 



The indices v and p indicate that the fat link Vi^ vp in direction fi is not decorated with staples extending 
in directions v and p. The projection to SU(3) can be defined in two different ways. W S SU(3) is the 
deterministic projection of A, 

W = Proj st/(3) A if Re Tr (WA f ) = Re Tr (UA*) , (2) 



*On leave from Institute for Theoretical Physics, Eotvos University, Budapest, Hungary. 

1 Note, that using a different projection one can calculate the fermion force 4 . This idea was applied for the fat link irrelevant 
clover (FLIC) action Recently a completely analytic smearing without any projection has also been introduced 6 . 
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whereas the probabilistic projection W x of A is chosen according to the probability distribution 



P(W) cx exp 



^ReTr^At) 



(3) 



In the second step the decorated links Vi tfl - V are constructed using the fat links Vi^ vp obtained in the first step 



as 



V, 



Proj s 



(7(3) 



In the final step the blocked links Vi „ are constructed as 



V, 



3) 



(4) 



(5) 



The smeared link obtained using the above construction containes thin links only from hypercubes attached to 
the original thin link. 

The HYP action is of the form 

S = S g (U) + S g (V) + S f (V), (6) 
where S g (U) is the plaquette gauge action 



depending on the thin links {£/} and S g (V) is the gauge action depending on the smeared links {V} [2J. The 
main role of S g (V) is to increase the acceptance rate in the accept- reject step. The simplest choice is the 
smeared plaquette S g (V) — — ^ ^ p Re Tr(V^), where 7 can be tuned to maximize the acceptance rate. The 
fermions are coupled to the smeared links, thus, the staggered fermionic matrix is of the form 



M(V)ij = 2m5 ii j + Vi,n (Vi^&ij-p, 



T/t X 

v i-p,,p, u ' l ,3+V 



(8) 



The matrix M(VyM(V) is block diagonal on even and odd lattice sites. Let denote the even block 

Q(V) := (M(V)<M(V)) eveneven . (9) 
Then the fermionic action Sf(V) describing rif flavors is of the form 



Sf(V) 



71/ 



Trlnfi(^). 



(10) 



Since the dependence of the smeared links {V} on the thin links {U} is nonlinear due to the projections 
to SU(3), the explicit form of the fermionic force, which is needed for molecular dynamics simulations, is very 
complicated. This fact makes the HMC and R algorithms virtually unusable. A two step algorithm, the partial 
global stochastic Metropolis update is used instead. In the first step a subset of the thin links {U} is updated 
such that the detailed balance condition with the thin link gauge action S g (U) is satisfied. This can be done 
using either heatbath or overrelaxation. In the second step the new smeared links {V} are computed and the 
newly obtained configuration is accepted with the probability 



F acc = min 1, exp [S g (V) + S g {V)] 



detO(F') 
detn(V) 



(11) 



Instead of calculating the ratio of the determinants a stochastic estimator is used. The ratio can be expressed 
as an expectation value 



det Q(V) 
dct Q(V) 



(exp(-r [SliVT^iV) - 1] £)> 



(12) 
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Only one random source £ is used on every gauge configuration pair {[/} and {£/'} to estimate the determinant 
ratio. The expectation value is taken together with the configuration ensemble average. Then the stochastic 
acceptance probability becomes 

P stoch = min{l,exp (-S g (V) + S g (V)) exp (-£* [Slfy'^QiV) - l] £)}. (13) 

If the stochastic estimator has large fluctuations then the acceptance rate can be very small even if the old 
and new fcrmionic determinants are almost the same. The standard deviation of the stochastic estimation 

dpt O' 

(exp (—ASf)) = = det(A)- 1 = (exp (-£* [A - 1] (14) 

can be written in the form 

a 2 = ( eX p(-2r[A-l]£)) e£ -(exp(-r[A-l]£))^ 

= dot (2,4 -I) -1 -det(A)- 2 , (15) 

where Q, = ft(V), fl' = Q(V') and A = fT^f! As suggested in QUI [TT], instead of fl and fi' we introduce the 
reduced matrices 

fi r = fie- 2 ^\ n' r = n'e- 2 ^ n '\ (16) 

where / is a polynomial chosen such that the function e 2 ^ x ^ jx is close to 1 in the interval where the eigenvalues 
of the matrix f2 can be found. Then the ratio of the determinants can be rewritten as 

= (exp(-r [^- 1 a.-l]0) r ? exp(2Tr [/(J)') - /(fi)]) . (17) 

Since the second factor can be evaluated exactly only the first factor has to be evaluated stochastically. Due 
to the special choice of /, A r = Sljr 1 ^!,- ss 1, so the fluctuations of the stochastic estimator are minimized, 
improving the acceptance rate. 

Equation (|15(l is valid for the standard deviation of the stochastic estimator only if the matrix 2A — 1 
is positive definite, that is, all eigenvalues of A are greater than 1/2. This is, however, very unlikely if the 
updating method changes a large number of links of the configuration. In order to avoid this problem the 
reduced fermionic determinant ratio can be written in the form 



det^r^det^V") n =(exp(-^e;[^ /,l -l]0]) - (18) 



where n is an arbitrary positive integer and £j are n independent random vectors. Then the standard deviation 
becomes 

a 2 = det 



(2A\l n -\\ -det (A r y 2 . (19) 



This is valid only if all the eigenvalues of A r are greater than 2~ n . If n is chosen large enough this condition 
can always be fulfilled. Since the determinant of a matrix product is independent of the order of the matrices, 
the nth root of A r can be written as 

A x J n = n'- 1 / 2n n 1 r / n n'~ 1 / 2n . (20) 

The factors can be approximated by polynomials as 



n i/n = nVn exp (_ 2/ (0)/n) = Q^ ) (0). (21) 



Here P^ 2 ™* 1 and are / and k order polynomials of the fermionic matrices f£ and f2', respectively. Then all 
the terms of the exponent of (|18|l can be written in the form 



r 



AV n - 1 



i = rp/ 2n) (o')Qi" ) (o)p/ 2n) (ooe - r^/ 2n) (wL n) (W/ 2n) (tt')e (22) 
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The polynomial orders I and k required for a reasonable approximation depend on the used quark mass. The 
polynomials should be optimized for the interval spanned by the eigenvalues of il. The smallest possible 
eigenvalue is 4m 2 , so if smaller and smaller quark masses are used then higher order polynomials are required. 
The polynomials have to be generated only once before each simulation using the method described in the 
following section. 



3 Generating the polynomials 



Our aim is to approximate the function / in the interval [a, b] using an nth order polynomial Pf (x). We 
choose a weight function w(x) and define the deviation of Pi n \x) from / using the distance in the Hilbert space 
C 2 w2 ([a,b]) as 



S n = 



f-P, 



(n) 



i-p ( ;\f-p ( ; ] 



f{x)~P { f n '(x) w{x)' 2 dx 



\f(x)\ 2 w(x) 2 dx 



(23) 



Here (,) w and || || denote the inner product 
and the norm 



f(x)*g(x)w(x) 2 dx 



Wf\L = y/(fJ)i 



(24) 



(25) 



in the Hilbert space L 2 w ^ ([a, b]), respectively. For the best choice of w(x) see Section I3~2l In order to minimize 
5 n we take a basis of orthogonal polynomials <E> M , 



where $ M {ji = 0, 1, 2, . . . ) is a /ith order polynomial with norm 



3m=II*X= f \®»{x)\ 2 w(x) 2 dx. 

J a 



(26) 



(27) 



This basis of polynomials is generated using the Gram-Schmidt orthogonalization process from the simple 
polynomials id (=1), id 1 (=x), id 2 (=x 2 ), . . . : 



$ 



id° = l 

u 



I 



|$o)($ol \ i i (^o^d 1 ) 
1 ° A ° iw id 1 = id 1 — i — - ^ $ 

ll*o|£ 

l*o)<#o| t 



'7o 

\ „ - (*o,id 2 ) ($i,id 2 ) 
1 id = id - — $q - — $i 



(lo 



'li 



1 - 1 id ' +1 = id ' +1 - 1 f <**,id'* 1 >.*,. 

V u=0 \\®A W ) v=0 qv 



Using this basis / can be written as 



where 



(28) 



(29) 



(30) 
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The nth order polynomial P^ for which S n is minimal can be obtained by taking only the first n terms of the 
sum in f2"3|) . 

In order to obtain the polynomials 3?^ a second order recursion formula |12l II 3j can be used instead of the 
numerically unstable Gram-Schmidt orthogonalization process (128(1 . The recursion goes as follows. The first 
two polynomials are obtained identically according to 1(28(1 

6 

w{x) 2 x dx 

$ Q (x) = l, $ l {x)=x- J -y b . (32) 

/ w(x) 2 dx 

J a 

Let 

p^:= / <S>^(x) 2 w(x) 2 xdx (33) 

J a 

and 

0„:=-^, 7.-1 :=— (34) 

9m 9m- i 

Then the rest of the polynomials can be obtained as 

^n+iix) = (x + p^Q^x) +7 M _i$ M _i(a;). (35) 



3.1 Numerical realization 

One proper method of proceeding with the algorithm is to calculate the integrals ((27(1 . ((30(1 and l(33|l exactly. 
This method is followed in Ref. |14|. These calculations require a precision of several hundreds or thousands 
of digits which can be carried out using multiprecision arithmetics libraries. Since the integrals are carried out 
exactly the indefinite integrals of the integrands have to be known. As a consequence only special types of 
functions / can be approximated and the weight function w also has to be carefully chosen. 

The method we use is to calculate the integrals numerically. The interval [a, b] is divided into N subinter- 
vals (N+l discretization points) logarithmically (see Section l3~2"]l . Since all the integrals are calculated using 
Simpson's rule, N has to be even. The values of the function /, polynomials <& M and the weight function w 2 
are calculated and stored only at the discretization points. First $o and <&i are determined with the integrals 
in 1(32(1 carried out numerically. In the [ith step and are calculated first using ((27(1 and 1(33(1 . then (3^ and 
7 M _i using 1(34(1 . Then is determined at each discretization point from $ M and $^-1 using ((35(1 . Then b^ 

and c M are calculated using ((30(1 . Finally, • $ M is added to the actual value of Pj 7 ^ . 

This method has three major advantages. Firstly, we have a second order recursion formula ((35(1 . Therefore, 
at each step only the last two orthogonal polynomials have to be stored in memory. That is, the memory 
required for the calculations depends only on the number of discretization points N but is independent of the 
order of approximation. Secondly, no multiprecision arithmetics is needed. All the calculations can be carried 
out using the built in 10 Byte floating point type. Finally, since the integrals are evaluated numerically, the 
indefinite integrals of the integrands are not needed. Therefore, there are no restrictions on the form of the 
function / and the weight function w. 



3.2 Stability, optimal weight function and discretization 

In order to describe the numerical stability of the recursion formula 1(35(1 and to find the optimal weight function 
w and the optimal type of discretization we need to refer to C 2 ([a, b]), the Hilbert space of [a, b] — ► C square- 
integrable functions with the inner product (/, g) — J f{x)*g(x)dx. If the weight function w is such that 

< K\ < w(x) < K2 Vs e [a, b] , (36) 

then the equivalence classes in C 2 ([a, b]) consist of the same functions as in ([a, b]) and C 2 ([a, b]) consist of 
the same equivalence classes as £^ 2 ([a, b]). That is, ([a, b]) and C 2 ([a, b]) are identical as linear spaces. 
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Both 



and 



{s„ nGN}u{c„ | n6N}u{co} 



form orthonormal bases in C 2 ([a, b]), where 



e n (x) 



exp 



lixni 
b — a 



a + b 



and 



s n (x) 

c n (x) 
Co (a;) 



./ 2 


V &- 


a 


./ 2 


V 6- 


a 


1 





2irn 
b — a 
27m 



\/b — a 

Using these basis vectors the linear map 



b — a 
e (x). 



a + b 
2 

a + b 



-,-2,-1,0,1,2, 

n = 1,2,3, ... 
n = 1,2,3,... 



F :=I + i- (s)J3 (l«n)(cn| - |cn)<«„|) = |e )(e | + 2 • (s)^ |e„)(e n | 



can be defined, which is bounded and self-adjoint in C 2 ([a, b\). Let 

A f (x) :=\(Ff)(x)\ 

and 

ip f (x) := arg ^(F/)(x) 

Then 



If / is real, then 



(Ff)(x)=Af{x)-exp{iip f {x)). 
f(x) = Re [(Ff) (x)^j = A f (x) ■ cos ( Vf {x)) . 



(37) 
(38) 

(39) 



(40) 



(41) 

(42) 
(43) 
(44) 
(45) 



That is, Af{x) and tpf(x) can be naturally identified as the amplitude and phase of the real function / at point 
x, respectively. If (pf is differentiable then we can define 



g f (x) := -tp'Ax). 

TT J 



(46) 



If / is a polynomial without multiple zeros (which is the case when dealing with orthogonal polynomials) then 
tpf is strictly increasing and Qf > 0. In this case if 



Qf(x) dx = k 



(47) 



for some a < X\ < X2 < b, then / has exactly k zeros in both [.ti,X2[ and jaji,^]. Thus, Qf can be identified 
as the density or root density of polynomial /. Graphically speaking, Af(x) describes the 'amplitude' of the 
polynomial / at point x and Qf{x) describes the 'rate at which the polynomial / oscillates' near x (Figure^). If 
/i and fi are polynomials without multiple zeros such that /i and fi have no common zeros, then the number 
of roots of /i • ji in every [xi , x^\ C [a, b] is equal to the sum of the number of roots of f\ and fi in that interval. 
Therefore, the root density of such a product is approximately the sum of the root density of the factors. 



G 




Figure 1: The amplitude (left), phase and density (right) of the 10th orthogonal polynomial Pio with respect 
to the weight function w — 1 on the interval [0, 4]. 



Let $ M and P M denote the /xth orthogonal polynomials in ([a, b]) and C 2 ([a, b]), respectively. If a = — 1 
and 6=1 then the polynomials are equal to the Legendre polynomials. Using the formulae for the asymptotic 
behaviour of the Legendre polynomials |15l §8.21] the formula 

A Plt (x) = const • ; 1 = + 0(/i" 3/2 ) (48) 

b — a\ 2 ( a + 6 X 2 



can be obtained for the amplitude Ap of the polynomials P M for large /is. 
If w satisfies condition (|36|) . then 

w{x)A^ ii (x) w const • A Pfi (x) (49) 
for large /is, that is, for large [is the amplitude of "I?^ can be well approximated by 

(x) w const • — !-y (x) , (50) 

where the constant is independent of x and is near 1 (Figure |2J) . Combining equations l|48|) and (|50|l the 
amplitude of the /ith orthogonal polynomial ^ can be approximately given by the formula 

A&Jx) w const • —5— ; =: (51) 

W(X) / / , x 2 / ,x 2 

v ' I ' b — a\ I a + cm 



For large /xs the root density of the polynomials can be approximated by 

1 

Q<s>^{x) m fiun- ; = n = =: fiu^- g(x), (52) 

b — a\ ( a + 6 X 
- \ x - 
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Figure 2: Left: A w , denned in equation l|51|) . as the amplitude of <I>60- Right: The relative deviation of the 
amplitude of <i>6o an d A w with a = 10~ 5 , 6 = 4 and w(x) = a; -1 / 4 . (The amplitude of $60 wa s calculated 
numerically using only the first 5000 terms of the sum in equation jjjl -) 




1 2 3 

x 

Figure 3: ^ • 0.9875 • g$ 60 (x) calculated numerically (solid line) such that only the first 5000 terms of the sum 
in equation (|41JI were taken into account and g(x) (dashed line) with a = 10~ 5 , b = 4 and w(x) — a; -1 / 4 . 
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0.01 



0.005 



-0.005 



-0.01 



i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — r 
F - $ 

1 w,60 *60 




J I I I I I I I I I I I I I I I I I 1_ 



2 
x 



Figure 4: Left: F w> eo (solid line) with uiqo = 0.9875 and $60 (dashed line). Right: The difference of F w eo and 
$ 60 . ( a = 1(T 5 , b = 4 and w(x) = x" 1 / 4 ) 



where u>a depends on the weight function w. In all cases w,, is close to 1 and lim u, L = 1 (Figure |3J). 

/J — >OG 

With 

<f(x) '.— 7r / glz) dz — 7r — arccos — ; (53) 

J a V b ~ a J 



and using equations fBTJ) and j^S} we can conclude at 



$ M (x) w A w (x) ■ cos <^ /iw 



7r — arccos 



2x — a — b 
b — a 



(54) 



an approximate formula for the /xth orthogonal polynomial with respect to the weight fuction w (Figure . 

We use Simpson's rule in order to calculate the integrals (|27|) . I|3t)|l and (|33|l numerically. The numerical 
integral of the function g taken over the discretization interval I = [z — ft, z + ft] is 



/.num. 



1 4 1 

-g(z - ft) + -g{z) + -g(z + ft) 



Assume that g is analytic in / with radius of convergence at z greater then ft, that is, 



g(z + u) = Y^ 



/")(z) 



— ft < u < ft. 



(55) 



(56) 



The exact integral and the numerical integral of g then becomes 



_ ~ g(">(z) ~ g(")(z) (-hr+l ~ gW\z) u2k+1 

g ^ n \ n+l n\ n + 1 ^ 

n=0 n=0 



and 



7,num. 



fe=0 



(2fc + l)! 



.9 = 



1™ z) 4^o (n) (z) l^a( n )(V) 



n=0 



n=0 



n=0 



ft = 



fE^W^) 



3 ^o ^ 



ft, 



(57) 



(58) 



9 



respectively. Then the error of integration over the interval /, that is, the difference of the exact and the 
numerical integrals becomes 



A/<? = 



/.num. 



E 

fc=2 



9 {2k) {z) (2 



[2k)\ V 3 2fc + ! 



2fe+l 



(59) 



If h is small enough and fi is large enough such that $^ oscillates much faster than w and / the integrand of 
(|30|1 can be approximated on / by 



F z (u) := B(z) cos (7r /i g(z) ■ u + Lp®^ (z)) — h < u < h, 



where B(z) = w(z) 2 f(z)A^ >ii (z). 



F^ k \u) = {-l) k ^ng{z)f k -F z {u), 



(60) 
(61) 



thus, the error of integrating F z numerically over / can be estimated as 



ArF, 



\F z (z)\ 



E 

k=2 



(-l) k (n^g(z)) 



2 A 



(2fe)! 



3 2fc + 1 



2A;+1 



<-/,|F z (z)|^ 

A:=2 



■2k ,2k 



(it \x g(z)) h 

m 



< ^h\B(z)\ 



cosh (tt /i g(z) h) — 1 — 



(TT^g(z)hy 



1 

36 



/i|5(z)| (tt^z)^ 4 • 1 + ((tt fj, g(z) h) 



(62) 



1 / (it n g(z) h) shows the number of discretization points between two consecutive zeros of $^ near z. The 
density of discretization has to be chosen such that the relative error 



5rF z 



AjF z 
\B(z)\-h 



(63) 



of the numerical integration is in the same order of magnitude in each discretization interval and is equal to the 
desired relative error of the numerical integration over [a, b]. 

Since cos 2 x = (1 + cos2a;)/2 and the integrands in (|2*7jl and ifHT^I contain the square of <5> p , they can be 
treated as a constant plus a cosine of double frequency in every discretization interval. The error of the numerical 
integration can be estimated similarly, leading us to the same conclusion. 

If / is not continuous on [a, b], then the convergence in (|29|l is not uniform but an almost everywhere 



convergence. Assume that / has an (x— a) a type (a > 0) singularity near a, that is, < 



lim f(x)(x — a) c 



< 



oo. Then the closer we get to a the slower the convergence becomes. Thus, by taking into account only a finite 
number of terms in (|29|) we can get a reasonable approximation for / only in [a + e, b] with a suitably chosen 
e. The procedure in this case should be as follows. We need to determine the size e of the neighborhood of a 
in which the approximation of / is not needed. Then we generate the orthogonal polynomials using H35J) in the 
interval [a' — a + e, b] and calculate the approximating polynomial Pj with the desired degree of approximation 
n. 

If / does not have singularities and is continuous on [a, b] then the weight function w can be chosen to be an 
arbitrary continuous function. If / has singularities, for example an (x — a)~ a type singularity near a, then the 
best choice for w is as follows, w should have an (x — a) a type behaviour near a and should be a smooth function 
otherwise. If / is such that \ f{x)\ > n > Vs e ]a,6] then the best choice is w(x) — 1/ |/(a;)|. Choosing w in 
such a way has the following advantages. 1. According to (|50[1 the amplitude of the polynomials <& M will gain 

(n) 

approximately the same type of singularity as / has, therefore, the relative deviation of / and is decreasing 
uniformly. 2. In the integrand of equation (|30J) the singularity of / is cancelled out by one of the two w's, while 
the other w deals with $ M according to (|50|l . 

In order to find the optimal type of discretization of interval [a, b) we need to consider both the singularities 
and the densities of the integrands (|27|) , l|33|l and l|30|) . Taking (|48|l and (|50|l into account we can conclude that 
|S(x)| of l(2T|) and l(3*3*|l (see (E2J)) has singularities near a and b of type \ j\Jx — a and 1/y/b — x, respectively, 



10 




-8 



10 



~i — i — i — i — jun 



-i 1 1 r 



N= 10000 



N=20000 

N=50000 



_l I I I I 1 I I I I L 



pi m 



•-•-.J 



2000 



4000 



6000 



cu 
I 



O 



-5 



10 



1 1 

- 


1 1 


1 i j i 1 1 


i-i--- 

r 




X^^^ i 








si 






-N=10000 - 






" N = 20000 - 






" N=50000 
, , , 1 


i i 1 i i i 


i , 



2000 



4000 
n 



6000 



Figure 5: The logarithm of the coefficients (left) and the deviations 



f-P 



(n) 



s 



(right) as the function of 



the order with number of discretization points N = 10000, N = 20000 and N = 50000. 



and (|30|l has singularities of type 1 / \[x — a and 1 / \fb — x. The densities of all three integrands are of type 
l/\/ {x — a)(b — x). As a consequence, the optimal discretization should contain the discretization points with 
density proportional to l/[(x — a)(b — x)]. This would require infinitely many discretization points, thus, we 
choose a small e and discretize the interval I' = [a' , b'] = [a + e, b — e] with the above density. If we use N 
discretization points, then the discretization interval [z — h z , z + h z ] at z has the length 



2h, = I 



(z — a) 
(b-z) 



b — a 
2e 

b — a 
2e 



4/N 



4/N 



a + b 



if a + e < z < 



if ^< z<6 _, 



The length of the longest and shortest intervals of this logarithmic discretization is 



4/N 



, 2e I 



1 



(64) 



and 



1 



respectively. 

In the usual fermionic calculations a = and the functions of the fermionic matrix that have to be evaluated 
have singularities of type x~ a (a > 0) at x — 0. We use the above polynomial expansion to approximate these 
singular functions, e should be chosen such that all eigenvalues of the fermionic matrix are greater then e. The 
smallest eigenvalue of the fermionic matrix is proportional to the square of the quark mass. Since our aim is 
to use quark masses as low as the physical u,d quark masses, which are approximately m u d = 0.002 in lattice 
mass units, e should be in the order of magnitude of 10~ 6 . In order to be able to well approximate the functions 
of the fermionic matrix so close to their singularities the required order of polynomial approximation is in the 
thousands. 

Using this logarithmic type of discretization the order n up to which the algorithm l|35|l is stable can be 



tested. The coefficients c„ and the deviations 



can be seen in Figure when the approximated 



function is chosen to be / = = 0, b = 4 and e = 10 6 . It can be verified that the algorithm is stable 

approximately up to the orders n = 1000, n = 2000 and n = 5500 if the numbers of discretization points are 
N = 10000, N = 20000 and N = 50000, respectively. Since no multiprecision arithmetics is required for the 
algorithm, its CPU time and memory requirements are considerably low. Generating the polynomials up to the 
order of n — 5500 in the case of TV = 50000 takes approx. 2.5 minutes of CPU time on a 1.6 GHz AMD Athlon 
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processor and requires approx. 5 MB of RAM. 



4 Conclusion 

We have introduced an alternate numerical method for generating the approximating polynomials used in 
fcrmionic calculations with smeared link actions. This algorithm was based on the idea of calculating all 
the integrals numerically and calculating and storing all the functions and polynomials only at finitely many 
discretization points. The advantages of this algorithm include memory usage independent of the order of 
approximation, unnecessarity of multiprecision arithmetics libraries and the absence of restrictions for the 
form of the approximated and the weight functions. We investigated the stability of the algorithm and based 
on the asymptotic behaviour of the orthogonal polynomial base appearing in the method we determined the 
optimal weight function and the optimal type of discretization. As a result the achievable order of polynomial 
approximation reached several thousands which is essential for fermionic calculations near the small physical 
quark masses. 
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